Status `~/Documents/github.com/ucla-biostat-216/2024fall/slides/14-app/Project.toml`
[1e616198] COSMO v0.8.9
[13f3f980] CairoMakie v0.12.16
[f65535da] Convex v0.16.3
[31c24e10] Distributions v0.25.113
[ced4e74d] DistributionsAD v0.6.57
[5789e2e9] FileIO v1.16.6
[f6369f11] ForwardDiff v0.10.38
[a2cc645c] GraphPlot v0.6.0
⌃ [86223c79] Graphs v1.9.0
[916415d5] Images v0.26.1
[ee78f7c6] Makie v0.21.16
[90014a1f] PDMats v0.11.31
[47aef6b3] SimpleWeightedGraphs v1.4.0
[0c5d862f] Symbolics v6.22.0
[e88e6eb3] Zygote v0.6.73
Info Packages marked with ⌃ have new versions available and may be upgradable.
1 Spectral clustering
Spectral clusting is one commonly used clustering approach in machine learning. ScikitLearn tutorial.
Given an undirected, weighted graph \(G=(V, E)\), we can define the weighted adjacency matrix\(\mathbf{W}=(w_{ij})_{i,j=1,...,n}\) with nonnegative weights \(w_{ij}\). \(w_{ij}=0\) means that the vertices \(v_i\) and \(v_j\) are not connected. The degree of a vertex \(v_i\) is \(d_i = \sum_{j=1}^n w_{ij}\).
The unnormalized graph Laplacian is \[
\mathbf{L} = \mathbf{D} - \mathbf{W},
\] where \(\mathbf{D} = \text{diag}(d_1,\ldots,d_n)\).
Exercise: For this graph, what is \(\mathbf{L}\)?
Properties of the unnormalized graph Laplacian:
For every vector \(\mathbf{f} \in \mathbb{R}^n\), we have \[
\mathbf{f}' \mathbf{L} \mathbf{f} = \frac 12 \sum_{i,j=1}^n w_{ij} (f_i - f_j)^2.
\]
Proof: BV exercises 3.21, 7.9.
\(\mathbf{L}\) is symmetric and positive semidefinite.
Proof: Part 1.
The smallest eigenvalue of \(\mathbf{L}\) is 0, the corresponding eigenvector is the constant one vector \(\mathbf{1}_n\).
Proof: Obvious.
(Number of connected components and the spectrum of \(\mathbf{L}\)) The multiplicity \(k\) of the eigenvalue 0 of \(\mathbf{L}\) equals the number of connected components \(A_1,\ldots,A_k\) in the graph. The eigenspace of eigenvalue 0 is spanned by the indicator vectors \(\mathbf{1}_{A_1}, \ldots, \mathbf{1}_{A_k}\) of those components.
Normalized graph Laplacians. There are two versions of normalized graph Laplacians in the literature: \[
\mathbf{L}_{\text{sym}} = \mathbf{D}^{-1/2} \mathbf{L} \mathbf{D}^{-1/2} = \mathbf{I} - \mathbf{D}^{-1/2} \mathbf{W} \mathbf{D}^{-1/2}
\] and \[
\mathbf{L}_{\text{rw}} = \mathbf{D}^{-1} \mathbf{L} = \mathbf{I} - \mathbf{D}^{-1} \mathbf{W}.
\]
Unnormalized spectral clustering algorithm.
Input: A similarity matrix \(\mathbf{S} \in \mathbb{R}^{n \times n}\) and number \(k\) of clusters to construct.
Construct a similarity graph. Let \(\mathbf{W}\) be its weighted adjacency matrix.
Compute the unnormalized Laplacian \(\mathbf{L}\).
Computer the first \(k\) eigenvectors \(\mathbf{u}_1, \ldots, \mathbf{u}_k\).
Let \(\mathbf{U} = (\mathbf{u}_1 \cdots \mathbf{u}_k) \in \mathbb{R}^{n \times k}\).
Treat rows of \(\mathbf{U}\) as data points in \(\mathbb{R}^k\) and cluster them using the \(k\)-means algorithm into clusters \(C_1,\ldots,C_k\). Output: Clusters \(A_1,\ldots,A_k\) according to \(C_1,\ldots,C_k\).
Normalized spectral clustering according to Shi and Malik (2000)
Input: A similarity matrix \(\mathbf{S} \in \mathbb{R}^{n \times n}\) and number \(k\) of clusters to construct.
Construct a similarity graph. Let \(\mathbf{W}\) be its weighted adjacency matrix.
Compute the unnormalized Laplacian \(\mathbf{L}\).
Computer the first \(k\) generalized eigenvectors \(\mathbf{u}_1, \ldots, \mathbf{u}_k\) of the generalized eigenproblem\(\mathbf{L} \mathbf{u} = \lambda \mathbf{D} \mathbf{U}\).
Let \(\mathbf{U} = (\mathbf{u}_1 \cdots \mathbf{u}_k) \in \mathbb{R}^{n \times k}\).
Treat rows of \(\mathbf{U}\) as data points in \(\mathbb{R}^k\) and cluster them using the \(k\)-means algorithm into clusters \(C_1,\ldots,C_k\). Output: Clusters \(A_1,\ldots,A_k\) according to \(C_1,\ldots,C_k\).
Remark: Generalized eigenproblem \(\mathbf{L} \mathbf{u} = \lambda \mathbf{D} \mathbf{U}\) is same as the eigenproblem \(\mathbf{L}_{\text{rw}} \mathbf{u} = \lambda \mathbf{u}\).
Normalized spectral clustering according to Ng, Jordan, and Weiss (2002)
Input: A similarity matrix \(\mathbf{S} \in \mathbb{R}^{n \times n}\) and number \(k\) of clusters to construct.
Construct a similarity graph. Let \(\mathbf{W}\) be its weighted adjacency matrix.
Compute the normalized Laplacian \(\mathbf{L}_{\text{sym}}\).
Computer the first \(k\) generalized eigenvectors \(\mathbf{u}_1, \ldots, \mathbf{u}_k\) of \(\mathbf{L}_{\text{sym}}\).
Let \(\mathbf{U} = (\mathbf{u}_1 \cdots \mathbf{u}_k) \in \mathbb{R}^{n \times k}\).
Form the matrix \(\mathbf{T} \in \mathbf{R}^{n \times k}\) from \(\mathbf{U}\) by normalizing the rows to norm 1.
Treat rows of \(\mathbf{T}\) as data points in \(\mathbb{R}^k\) and cluster them using the \(k\)-means algorithm into clusters \(C_1,\ldots,C_k\). Output: Clusters \(A_1,\ldots,A_k\) according to \(C_1,\ldots,C_k\).
2 Matrix completion
Snapshot of the kind of data collected by Netflix. Only 100,480,507 ratings (1.2% entries of the 480K-by-18K matrix) are observed
Netflix challenge: impute the unobserved ratings for personalized recommendation. http://en.wikipedia.org/wiki/Netflix_Prize
Matrix completion problem. Observe a very sparse matrix \(\mathbf{Y} = (y_{ij})\). Want to impute all the missing entries. It is possible only when the matrix is structured, e.g., of low rank.
Example: Load the 128×128 Lena picture with missing pixels.
We fill out the missing pixels uisng a matrix completion technique developed by Candes and Tao \[
\text{minimize } \|\mathbf{X}\|_*
\]\[
\text{subject to } x_{ij} = y_{ij} \text{ for all observed entries } (i, j).
\] Here \(\|\mathbf{M}\|_* = \sum_i \sigma_i(\mathbf{M})\) is the nuclear norm. In words we seek the matrix with minimal nuclear norm that agrees with the observed entries. This is a semidefinite programming (SDP) problem readily solved by modern convex optimization software.
We use the convex optimizaion package COSMO.jl to solve for this semi-definite program.
# # Use COSMO solverusingConvex, COSMOsolver = COSMO.Optimizer()# Linear indices of obs. entriesobsidx =findall(Y[:] .≠0.0)# Create optimization variablesX =Variable(size(Y))# Set up optmization problemproblem =minimize(nuclearnorm(X))problem.constraints += X[obsidx] == Y[obsidx]# Solve the problem by calling solve@timesolve!(problem, solver) # fast
------------------------------------------------------------------
COSMO v0.8.8 - A Quadratic Objective Conic Solver
Michael Garstka
University of Oxford, 2017 - 2022
------------------------------------------------------------------
Problem: x ∈ R^{49153},
constraints: A ∈ R^{73665x49153} (73793 nnz),
matrix size to factor: 122818x122818,
Floating-point precision: Float64
Sets: ZeroSet of dim: 40769
DensePsdConeTriangle of dim: 32896 (256x256)
Settings: ϵ_abs = 1.0e-05, ϵ_rel = 1.0e-05,
ϵ_prim_inf = 1.0e-04, ϵ_dual_inf = 1.0e-04,
ρ = 0.1, σ = 1e-06, α = 1.6,
max_iter = 5000,
scaling iter = 10 (on),
check termination every 25 iter,
check infeasibility every 40 iter,
KKT system solver: QDLDL
Acc: Anderson Type2{QRDecomp},
Memory size = 15, RestartedMemory,
Safeguarded: true, tol: 2.0
Setup Time: 35.32ms
Iter: Objective: Primal Res: Dual Res: Rho:
1 -1.4426e+03 1.1678e+01 5.9856e-01 1.0000e-01
25 1.4495e+02 5.5360e-02 1.1033e-03 1.0000e-01
50 1.4754e+02 1.2369e-02 1.6744e-03 7.4179e-01
75 1.4797e+02 4.9490e-04 5.5696e-05 7.4179e-01
100 1.4797e+02 1.4115e-05 1.3438e-06 7.4179e-01
------------------------------------------------------------------
>>> Results
Status: Solved
Iterations: 100
Optimal objective: 148
Runtime: 0.857s (856.75ms)
1.242270 seconds (1.73 M allocations: 675.980 MiB, 5.35% gc time, 26.32% compilation time: 97% of which was recompilation)
colorview(Gray, X.value)
3 Compressed sensing
Compressed sensingCandes and Tao (2006) and Donoho (2006) tries to address a fundamental question: how to compress and transmit a complex signal (e.g., musical clips, mega-pixel images), which can be decoded to recover the original signal?
Suppose a signal \(\mathbf{x} \in \mathbb{R}^n\) is sparse with \(s\) non-zeros. We under-sample the signal by multiplying a (flat) measurement matrix \(\mathbf{y} = \mathbf{A} \mathbf{x}\), where \(\mathbf{A} \in \mathbb{R}^{m\times n}\) has iid normal entries. Candes, Romberg and Tao (2006) show that the solution to \[
\begin{eqnarray*}
&\text{minimize}& \|\mathbf{x}\|_1 \\
&\text{subject to}& \mathbf{A} \mathbf{x} = \mathbf{y}
\end{eqnarray*}
\] exactly recovers the true signal under certain conditions on \(\mathbf{A}\) when \(n \gg s\) and \(m \approx s \ln(n/s)\). Why sparsity is a reasonable assumption? Virtually all real-world images have low information content.
Generate a sparse signal and sub-sampling:
usingCairoMakie, Makie, Random# random seedRandom.seed!(123)# Size of signaln =1024# Sparsity (# nonzeros) in the signals =10# Number of samples (undersample by a factor of 8) m =128# Generate and display the signalx0 =zeros(n)x0[rand(1:n, s)] =randn(s)# Generate the random sampling matrixA =randn(m, n) / m# Subsample by multiplexingy = A * x0# plot the true signalf =Figure()Makie.Axis( f[1, 1], title ="True Signal x0", xlabel ="x", ylabel ="y")lines!(1:n, x0)f
Solve the linear programming problem.
# Use COSMO solversolver = COSMO.Optimizer()# MOI.set(solver, MOI.RawOptimizerAttribute("max_iter"), 5000)# Set up optimizaiton problemx =Variable(n)problem =minimize(norm(x, 1))problem.constraints += A * x == y# Solve the problem@timesolve!(problem, solver)
# Display the solutionf =Figure()Makie.Axis( f[1, 1], title ="Reconstructed signal overlayed with x0", xlabel ="x", ylabel ="y")scatter!(1:n, x0, label ="truth")lines!(1:n, vec(x.value), label ="recovery")axislegend(position =:lt)f
4 Automatic differentiation (Auto-Diff)
Last week we scratched the surface of matrix/vector calculus and chain rules. Recent surge of machine learning sparked rapid advancement of automatic differentiation, which applies chain rule to the computer code to obtain exact gradients (up to machine precision).
4.1 Example: multivariate normal
Define log-likelihood function.
usingDistributions, DistributionsAD, LinearAlgebra, PDMats, Random, ZygoteRandom.seed!(216)n, p =100, 3Y =randn(p, n) # each column of Y is a sample
Surprise! Gradients of the covariance matrix do not match. Close inspection reveals that Julia calculates the gradient with respect to the upper triangular part of the covariance matrix \(\text{vech}(\Omega)\)
